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Abstract 

In recent works several authors have proposed the use of precise bound- 
ary conditions (BCs) for blurring models and they proved that the result- 
ing choice (Neumann or reflective, anti-reflective) leads to fast algorithms 
both for deblurring and for detecting the regularization parameters in 
presence of noise. When considering a symmetric point spread function, 
the crucial fact is that such BCs are related to fast trigonometric trans- 
forms. 

In this paper we combine the use of precise BCs with the Total Vari- 
ation (TV) approach in order to preserve the jumps of the given signal 
(edges of the given image) as much as possible. We consider a classic fixed 
point method with a preconditioned Krylov method (usually the conju- 
gate gradient method) for the inner iteration. Based on fast trigonometric 
transforms, we propose some preconditioning strategies which are suitable 
for reflective and anti-reflective BCs. A theoretical analysis motivates the 
choice of our preconditioners and an extensive numerical experimenta- 
tion is reported and critically discussed. The latter shows that the TV 
regularization with anti-reflective BCs implies not only a reduced analyt- 
ical error, but also a lower computational cost of the whole restoration 
procedure over the other BCs. 
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1 Introduction 

We are concerned with specific linear algebra/matrix theory aspects of the vast 
field of inverse problems [T71 [TB] which model the blurring of signals and images 
{2D or dD with d > 3). Here the goal is to reconstruct the real object from 
its blurred and noisy version and this goal is a classical one in astronomical 
imaging, medical imaging, geosciences, etc. [5]. 

The blurring model is assumed to be space- invariant, i.e., the point spread 
function (PSF) is represented by a specific bivariate function h{x — y) {x,y G il) 
for some univariate function h{-) [TB]. According to the linear models described 
in the literature |16) . the observed signal or image v and the original signal or 
image u are described by the relation 

v{x) = 'Hu{x) + ri{x) := / h{x — s)u{s)ds + ri{x), x £ (1) 
Jn 

where the kernel h is the PSF and r] denotes the noise. The problem ([T]) 
is ill-posed since the operator T-L is compact [TB]. Therefore, the approxima- 
tion/discretization matrix of H is usually increasingly ill-conditioned when the 
number n of pixels becomes large. In addition, the size of the subspace associ- 
ated with small eigenvalues, which substantially intersects the high frequencies, 
is large and proportional to the size of the matrix. Thus, we cannot directly 
solve Hu = d, since the small perturbations, represented by the noise 77 with 
important high frequency components due to its probabilistic nature, would be 
amplified unacceptably. 

To remedy to the latter essential ill-conditioning of problem ([I}, one may 
employ regularization methods. The Total Variation (TV) regularization ap- 
proach is a good choice for restoring edges of the original signals [22^ . Rudin, 
Osher, and Fatemi gave the total variation functional in the form 

Jtv{u) := / \Vu\dx, (2) 
Jn 

where | • | denotes the Euclidean norm. We note that the Euclidean norm | • | is 
not differentiable at zero. To avoid the non-differentiability, Acar and Vogel [5] 
considered the following minimization 

niin - v\\mn) + y/WW^ dx^ > (3) 

where a, /? are positive parameters. Notice that the penalty term -^/jVup -I- j3'^ dx 
converges to Jtv{u) as /3 0. In other words, the latter is a differentiable reg- 
ularized version of Jtv (u) ■ The corresponding Euler-Lagrange equation for ([3]) 
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is given by 

(4) 

fi = 0, xedn, 

on ' ' 

where * denotes the adjoint operator and Cu(y) —V • I , ^ Vy ) is the 

differential operator appearing in ^ and conies from the regularized penalty 
term. Vogel and Oman (29j proposed a lagged diffusivity fixed point (FP) 
iteration for solving Q. More precisely, given the initial guess u'^, the new 
iterate is obtained by m*^ thanks to the equation 

A,^ku''+^ = {H*H + aL{u''))u''+^ = H*v, fc = 0,l,..., (5) 

where H and L{u'^)u'^^^ denote the discretization/approximation matrix of H 
and £„fe (u*'"'"^), respectively. We use the compound mid point rule and stan- 
dard centered finite differences of precision order two for the finite dimensional 
approximation of H and Cu{-), respectively. Therefore, at each FP iteration, we 
may use the preconditioned conjugate gradient (PCG) method [EJ Algorithm 
10.3.1] for solving the linear system ([5]). 

In [5] the authors proposed a cosine preconditioner when iJ is a Toeplitz 
matrix, i.e., in the case of zero-Dirichlet boundary conditions (BCs). However, 
the choice of such BCs induces remarkable pathologies in the quality of the 
restored images, which should be avoided or at least minimized. In reality, using 
classical BCs such as periodic or zero-Dirichlet may imply, when the background 
is not uniformly black, disturbing Gibbs phenomena called ringing effects [201 

[21 EH]- 

The novelty of this paper is represented by the choice of appropriate BCs, in 
order to reduce the ringing effects, and in the related matrix/numerical analysis. 
The latter will affect the algebraic expression of H, while for L{u^) the choice 
of the BCs in the Euler-Lagrange equation for ^ seems to impose Neumann 
BCs. In this context, the idea is to combine the application of anti-reflective 
BCs, already studied for their precision with plain regularization methods like 
Tikhonov and Landweber [Ml [lH [Ml HH H [10] , with the more sophisticate TV 
regularization. In other words, the first aim consists in checking how to reduce 
the ringing effects and the over-smoothing of the edges simultaneously. Next, 
we want to study the use of preconditioners based on innovative fast trans- 
forms in the setting of Krylov methods when the real problem is modeled by a 
symmetric PSF. The final goal is to combine the precision of the reconstruction 
with highly efficient numerical procedures. We study some preconditioning tech- 
niques and give the theoretical explanations of different proposals. An effective 
preconditioner for the refiective BCs is inspired by the work in [5], while, for the 
anti-reflective BCs, we propose a new sine preconditioner for the linear system 
([5]) and explore the re-blurring approach introduced in pjj. Numerical results 
confirm the effectiveness of the proposed preconditioners and the superiority of 
the antirefiective BCs with respect to the refiective BCs. Indeed, antirefiective 
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BCs not only provide better restorations as expected (see [Ml E])) but also 
require a lower computational cost. The latter is a consequence of the fact that 
a more precise model requires lesser regularization, i.e. a smaller a, and the 
convergence of our approach is fast for small a. 

The paper is organized into seven more sections. In Section [5] we consider 
reflective and anti-reflective BCs. Sections[3]is devoted to define optimal precon- 
ditioners for signal and image deblurring for the reflective BCs and anti-reflective 
BCs. In Section |3] some spectral features of the proposed preconditioners are 
discussed. Section [5] is concerned with the numerical tests for checking the real 
efficiency of the considered preconditioners and the quality of the restorations. 
Finally, in Section [5] we draw conclusions. 



2 Boundary Conditions 

We start by introducing the one-dimensional deblurring problem. Consider the 
original signal u = (. . . , w-m+i, . . . , uq, ui, . . . , w„, w„+i, . . . , w„+m,, . . .)'^ and 
the normalized blurring PSF given by 



ho,---, 0,0,... y 



(6) 



with the usual normalization, i.e., X]j=-m ^ ^ which preserves the global 
intensity and therefore represents an average. The blurred signal v is the con- 
volution of h and u and consequently Vi — X]'i^-oo ^jUi-j is such that 



ho 



ho 



ho 



h-rn / 



U-m+1 
U-m+2 

Uo 

u 

Un+1 



(7) 



The deblurring problem is to recover the vector u — (ui , . . . ,Un) given the blur- 
ring function h and a blurred signal v = (vi, . . . , of finite length, Thus the 
blurred signal v is determined not only by u, but also by {u-m+i, - - - , uo)'^ and 
(u„+i, . . . ,Un+m)'^ and the linear system ^ is underdetermined. To overcome 
this, we make certain assumptions, that is the BCs on the unknown bound- 
ary data U-m+i, - ■ - ,uo and Un+i, - - - , Un+m in such a way that the number of 
unknowns equals the number of equations. 

For the zero-Dirichlet BCs, we assume that the data outside u are zero, 
i.e., we set ui-j = Un+j = for j ~ 1, . . . ,m. Then, ([7]) becomes Au — v, 
where A is Toeplitz. For the periodic BCs, we assume that the signal u is 



extended by periodicity. More precisely, we set = u„_ 



-i+i 



and 



''71+ j 
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for j = 1, . . . ,m. It follows that ([7]) becomes Au = v, where A is circulant 
and hence it can be diagonalized by the Discrete Fourier Transform (DFT) (see 

m)- 

For the Neumann or reflective BCs, we assume that the data outside u are 
a reflection of the data inside u (refer to [ID])- More precisely, we set ui-j = Uj 
and Un+j — it„+i_j for all j = 1, . . . , m in ([7]). Thus ([7]) becomes Au = v, where 
A is neither Toeplitz nor circulant but a special n-by-n Toeplitz plus Hankel 
matrix which is diagonalized by the discrete cosine transform provided that the 
blurring function h is symmetric, i.e., hj ~ h^j for all j in It follows that 
the above system can be solved by using three fast cosine transforms (FCTs) in 
O(nlogn) operations [50]. This approach is computationally interesting since 
the FCT requires only real operations and is about twice as fast as the FFT 
and this is true in two dimensions as well. We note that the reflection ensures 
the continuity of the signal and the error is linear in the discretization step 
(the latter can be easily seen by applying a Taylor expansion [10] )■ Therefore, 
we usually observe a reduction of the boundary artifacts with respect to zero- 
Dirichlet and periodic BCs. 

For the anti-reflective BCs, we assume that the data outside u are an anti- 
reflection of the data inside u. More precisely, if x is a point outside the domain 
and X* is the closest boundary point, then we have x = x* —Sx and the quantity 
u{x) is approximated by u{x*) — {u{x* + dx) — u{x*)). Consequently, we set 



ui-j = ui — (mj+i — til) = 2ui — Mj+i, for all j = \, . . . ,m, 

Un+j = Un - {Un-j - Un) = 2u„ - Un-j, for aU j = 1 , . . . , TO 



(8) 



in ^ . By a Taylor expansion, the anti- reflection ensures a continuity of 
the signal and the error is quadratic in the discretization step p¥ . Usually, the 
boundary artifacts are reduced also with respect to reflective BCs. 

Imposing the anti-reflection ([8]), the linear system ([7]) becomes Au — v, 
where ^ is a Toeplitz plus Hankel plus a rank-2 correction matrix, where the 
correction is placed at the first and the last column. Furthermore, in [3] the 
authors proved that if h is symmetric then A = TnAT~^ where 



1 

P Sn-2 JP 
1 



1 

-Sn-2P Sn-2 -Sn~2Jp 

1 



where J is the flip matrix, Sn-2 is the sine transform matrix of order n — 2, and 
where pj — 1 — j/(n — 1) so that the first column vector is exactly the sampling 
of the function 1 — 2: on the grid j /{n — 1) for j = 0, . . . , n — 1. Finally, A is a 
diagonal matrix given by suitable samplings of the function 



which is the symbol generated by the PSF. That is 
A = diagj^^i^ „ {h{yj)) , Sn, = 



2 I . I ]™ 

sm 



TO + 1 



TO + 1 



(9) 



(10) 
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where 

{j - l)n 

Uj = —, for j = 1, . . . , n - 1, and j/„ = 0. (11) 

n — 1 

As a consequence, a generic system Au = v can be solved within 0{n log n) 
real operations by resorting to the application of three fast sine transforms 
(FSTs) (refer to [23]), where each FST is computationally as cheap as a generic 
FCT. There is a suggestive functional interpretation of the transform T„. The 
transform associated with periodic BCs matrices is the Fourier transform: its j- 
th column vector, up to a normalizing scalar factor, can be viewed as a sampling, 
over a suitable uniform gridding of [0, 27r], of the frequency function exp(— ijy). 
Analogously, when imposing reflective BCs with a strongly symmetric PSF, the 
transform of the related reflective BCs matrices is the cosine transform: its j-th. 
column vector, up to a normalizing scalar factor, can be viewed as a sampling, 
over a suitable uniform gridding of [0,7r], of the frequency function cos{jy). 
Here the imposition of the anti-reflective BCs, by operating a central symmetry 
with respect to any point of frontier, can be functionally interpreted as a linear 
combination of sine functions and of linear polynomials (whose use is exactly 
required for imposing the continuity at the borders) . This intuition becomes 
evident in the expression of Tn ■ Indeed 

T„ = ( 1 - - , sm{y) , • • • , sin((n - 2)y) , ^ ) • diag ll,J In-2 , 1 

(12) 

where y is defined (|TTt and it is a suitable gridding of [0,7r]. 

Now, we introduce the two-dimensional case. For the Neumann or reflective 
BCs, the blurring matrix is a block Toeplitz-plus-Hankel matrix with Toeplitz- 
plus-Hankel blocks and can be diagonalized by the two-dimensional FCTs (which 
are tensor products of one-dimensional FCTs) in 0{n^ log n) operations provided 
that h is quadrantally symmetric, i.e., hi j = h-i j = hi -j = h-i (refer to 

m)- 

For the anti-reflective BCs, we assume that the data outside u are an anti- 
reflection of the data inside u, i.e., a point outside the domain is anti-reflective to 
the closest boundary point first in one direction and then in the other direction. 
In particular, we set 

= 2-ui,0 — = 2u„,0 — Mn-j,0, for 1 < j < m, 1 < (/) < n, 

— 2'«i/>,l — MiAj + l, U^,n+j — '2.U^,n — U^,n~], ioi 1 < j < Ul, 1 < 1p < 71. 

When both indices lie outside the range {1, . . . , n} (this happens close to the 4 
corners of the given image), we set 



Ml 




= 4ui4 


— 2uij-_ 




- 2Mi+i,i 4 


- Ui+ 




Ul- 




= 4ui,„ 


- 2mi,„. 


-j 


— 2uj+i^„ - 


t- Uti 


-l.n-j, 


Un- 


H.l-j 


= 4u„4 




n 


2ii^_2,i 


f- Un- 






-i.n+i 






-j 


2i/j2_.j_^ 


f Un 





for I < i, j < m. If the blurring function (PSF) h is quadrantally symmetric, 
then the blurring matrix is a block Tocplitz-plus-IIankel-plus-2-rank-correction 
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matrix with Toeplitz-plus-Hankel-plus-2-rank-correction blocks and can be di- 
agonalized by the two-dimensional anti-reflective transforms (which are tensor 
products of one-dimensional anti- reflective transforms r„) in 0(n^ log n) real 
operations (see for instance |3]). In the following we will assume a symmetric 
(quadrantally symmetric in 2D) PSF since reflective and anti-reflective BCs can 
be diagonalized by fast transforms only in such case. However, in the nonsym- 
metric case, even if the blurring matrix can not be diagonalized by fast trans- 
forms, the matrix- vector product can be done again in 0(n^ log n) by FFTs. 
Moreover, many practical blur have the symmetry like the celebrated Gaussian 
blur widely used in several contexts. 

3 Optimal Preconditioners with different Bound- 
ary Conditions 

The optimal preconditioner for a matrix A aims to find an approximation which 
minimizes A|| over all i? in a set of matrices for the matrix Frobenius norm 
II • II the typical set of matrices is formed by considering an algebra of matrices 
which are simultaneously diagonalized by a given unitary transform. The main 
novelty in our context is represented by the fact that the anti-reflective matrices 
with symmetric PSFs form a commutative algebra associated to a non-unitary 
transform. The latter poses nontrivial difficulties that are treated in the sequel 
of the paper. The optimal circulant preconditioner was originally given in [9]. 
The optimal sine transform preconditioner was presented in [7]. The optimal 
cosine transform preconditioner was provided in [^. In this section, we con- 
struct the optimal reflective BCs preconditioner and the optimal anti-reflective 
BCs preconditioner for ([S|) and the optimal reblurring preconditioner for the 
reblurring equation ^U\\ below instead of ([S]). Some of these preconditioning 
techniques are inspired from the idea proposed in [6l |8] for zero-Dirichlet BCs. 

3.1 One-dimensional Problems 

For the one-dimensional problems, we assume a symmetric and normalized PSF. 
Suppose that we impose the reflective BCs on H and the zero Neumann BCs 
on L{u^). In this case, we propose the following reflective BCs preconditioners 
for ([S]). Let Cn be the n dimensional discrete cosine transform with entries 

[C„],j = ^ — cos — y i,j = l,...,n, 

where 5ij is the Kronecker delta. The matrix C„ is orthogonal, i.e., CnC^ = I- 
Moreover, for any n-vector w, the matrix-vector product C'nW can be computed 
within 0(n log n) real operations by the FCT. Define 

C — {C^ACn : A is a real n-by-n diagonal matrix}. 
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For an n-hy-n matrix A, the optimal cosine transform preconditioner is 

c{A) = argmin \\B — A\\p. 
Bee 

The operator c(-) is hnear, preserves positive definiteness, and compresses every 
unitarily invariant norm. For the specific cosine algebra and a general conver- 
gence theory based on the Korovkin theorems, we refer to [H [23]. As in [8] 
for Dirichlet BCs, the optimal cosine transform preconditioner (i.e., the optimal 
reflective BCs preconditioner) for ([5]) can be defined as 

R = H*H + ac{L{u^)). (13) 

We note that R = c{A^k) since H* = H Q C. Spectral properties of the 
preconditioner will be discussed in Sectional Here, we only note that c{L{u^)) 
is not an optimal preconditioner for L{u^) if the coefficient (jVitp + /3^)^^^^ has 
large variation. In such case a diagonal scaling is necessary to obtain an effective 
preconditioner like diag(i(u'^))^/^c(L(it'^))diag(i(u'^))^/^, where diag(i(u'^)) is 
the diagonal matrix whose diagonal entries are the same as that of L{u^) pS] . 
We note that the coefficient matrix in ([5|) is the sum of two operators. To avoid 
the possibly large fluctuation in the coefficient of the operator in (O, we define 
a reflective BCs preconditioner for ([5]) by Dn = D^RD^^ where R is given in 
(pl) and 

D = / + adiag(L(u'=)). (14) 

A further possibility is to employ a diagonal scaling for ([5]). As in [8], we concern 
the scaled equation 

= (^H*H + aL{u'')^ {t'^+i = H*v, (15) 

where H ^ HD'^/^, L{u'') = D-^/^L{u'')D-^/^, and u'' = D^/^u''. Then, 
we propose the following reflective BCs preconditioner for (fTSj) Rd = H*H + 
ac{L{u^))^ where H — Hc{D~^^'^). If Ajj, A^i and A^ denote the eigenvalue 
matrices of c{D^^/'^) and c{L{u'^)), respectively, then Ro can be written as 

Rd = C^iAHAHA*D^D+aAi)Cn. 

Next, we construct the anti-reflective BCs preconditioners for ([5]) under the anti- 
reflective BCs for H and the Neumann BCs or anti-reflective BCs for L{u^). 
Let Sn be the n dimensional discrete sine transform of type I with entries as 
in PU)). Then, Sn is orthogonal and symmetric, i.e., = Sn and Sn = I- 
Moreover, for any n dimensional vector w, the matrix-vector product SnW can 
be computed in O(nlogn) real operations by the FST. Define r — {S'„AS'„ : 
A is a real diagonal matrix of order n}. Let ct{z) := {z2, ■ ■ ■ , Zn, 0)'^ with z = 
(zi, . . . , 2„)"^. Let T{z) be the n-hy-n symmetric Toeplitz matrix whose first 
column is z and 'H{z,Jz) be the n-by-n Hankel matrix whose first and last 
column are z and Jz, respectively. It was shown that for any B € t, there 
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exists z = (zi, . . . , z,,Y e such that [7] B = T(z) ~ H{cr^{z), Jcr'^{z)). For 
an n-hy-n matrix A, the optimal sine preconditioncr is 

s{A) = argmin \\B - A\\f. (16) 

BET 

The construction of s{A) requires only 0{n^) operation for a general matrix A 
and 0(n) operation for a banded matrix A. Furthermore, s(-) is linear, preserves 
positive definiteness, and compresses any unitarily invariant norm (see [71 123] V 
Now, we define an optimal sine transform based preconditioner (i.e., the 
so-called anti-reflective BCs preconditioner) for ([SJ by 

M = s{Hy s{H) + a s{L{u'')) (17) 

in the sense that, for any n-hy-n matrix A, s{A) is given by 

s{A) = aigmmWB - A\\f, (18) 

where = S'„AS'„ : A is a real diagonal matrix of order n and Sn :— diag(l, Sn-2, 1) 



Proposition 1 Given an n-by-n matrix A, we have 

KA)-- 



A(l,l) 

s(A(2 : n- 1,2 : n- 1)) 
A{n,n) 



where s(-) and s(-) are defined in (jl8p and (|16p . respectively, and A{2 : n— 1,2 : 
n — 1) is the suhmatrix of A corresponding to rows indexed from 2 to n—1 and 
columns from 2 to n — \ . 

Proof: By unitary invariance of the Frobenius norm [Sn is unitary) we find 
11^ — Sn^SnWr — llS'ri^'S'n — ^\\f, whcrc A is a diagonal matrix. To conclude 
the proof, it is enough to observe that 



di&g{SnASn) = 
and that s{A) = S'„diag(S'„74S'„)S'„ 



A{l,l) 

diag(S'„_2^(2 : 71- 1,2 : 71- 1)S'„_2) 
A{n,n) 



□ 



To reduce the potential fluctuations in the coefficient of the elliptic operator 
in ([5]), based on the diagonal scaling D in (fT4|) . we define a scaled anti-reflective 
BCs preconditioner for ^ by Dm — D^MD^^ where M is deflned in ([TT)). 
Similarly, for the scaled equation in the form of ()15p , we give the anti-reflective 
BCs preconditioner with diagonal scaling as follows Mo = H*H -\- as{L(u^))^ 
where H — s{H)s{D^^/^). If Ah, Ad and A^ denote the eigenvalue matrices 
of s{H), s{D^^/^) and J(L(m'^)), respectively, then Md can be written as 
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Finally, we consider the reblurring method with some new reblurring precon- 
ditioners under the anti-reflective BCs for H and the Neumann BCs for L{u^). 
In Section [21 we have observed that anti-reflective BCs matrices can be diago- 
nalized by the anti-reflective transform T„. Hence, it is possible to define the 
anti-reflective algebra 

ATZ = {T„Ar,7^ : A is a real diagonal matrix of order n\. (19) 

Unfortunately, H g ATZ but H* ^ AR. However, in [11], it was proposed to 
use a reblurring approach, i.e., to replace H* with H' , where H' is the matrix 
obtained by imposing anti-reflective BCs to the PSF rotated by 180 degrees. 
Since the PSF is assumed to be symmetric, H' = H [T2]- Therefore, instead of 
([S]), one may solve the following equation [TT] 

A'^^v!'+^ = {H'H + aL{u^))u^+^ ^H'v; /c = 0,l,... (20) 

by the PBiCGstab method |28| since A'^^. is not symmetric. In this case, a 
reblurring preconditioner for (j20p is given by 

P = H' H + a AR{L{u^)). 

Here, for any n-by-n matrix A, AR{A) is defined by 

zi + '^Yllhk •••0 



AR[A) 



Zn-2 

2zn-2 s{A{2:n-l,2:n-l)) Zn-3 + 2zn-2 



Zn~2 '■ 

^2 + 2 J^lZl Zk 

•••0 zi+2J2lZlzk 

where z = (zi, Z2, . . . , Zn-2)'^ is such that s{A{2 : n — 1, 2 : n — 1)) = T{z) — 
H{a'^{z), Ja'^{z)). We only need form s{A{2 : n — 1, 2 : n — 1)) for computing 
AR{A). 

We note that AR{A) belongs to the algebra ATZ defined in ([T9|. where A 
is defined as in (jlOp . Therefore, a linear system Au = v can be solved within 
0(n log n) real operations by using three FSTs. 

Remark 2 In general, AR{A) 7^ argmin^g^^^ — Ajji?. Moreover, we can not 
construct argmin^g^^^ \\B — A\\f in only 0(n^) operations by using the similar 
technique for computing s{A) for a general matrix A in Notice that Tn in 
dni) is Tn = Sn{I + U) and T'^ = {I - U)Sn, where 



U 



As in [19], we can compute the eigenvalue of ATl{A) by using the diagonal 
entries of ^ = SnASn, However, it requires O(n^logn) operations to calculate 
the diagonal entries of ^ . 
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To reduce the potential fluctuations in the coefflcient of the elhptic operator 
in (|20l) , we define a diagonally scaled reblurring preconditioner for ([20]) as follows 
Dp = D^/'^PD^/'^, where D is defined as the same form in ([Ti)) . For the scaled 
system 

A'^uu^+^ ^ {h'H + aL{u'')^ u'^+i = H'v, (21) 

the reblurring preconditioned is given by 

Pd = AR{D~^/^)H'H AR{D--^/'^) + aAR{L{u^)). 

If Kh, Ac, and denote the eigenvalue matrices of i7, AR{D^^/'^)^ and 
AR{L{u^y) , respectively, then the preconditioner Pd can be written as 

Pd = T,,{K*hAh^*d^d + aAL)T~\ 

A further possibility is the use of anti-reflective BCs for L{u''). This implies 
that the coefficient matrix in the linear equation ([^(7)) is closer to the precondi- 
tioner. Consequently, a faster convergence and a lower global cost have to be 
expected. The latter choice is in fact considered in the numerics. 

We comment on the cost of constructing Xjj, X E {R,M,P} and of each 
PCG/ PBiCGstab iteration. We note that L{u'') is a banded matrix. There- 
fore, computing c{L(u'^)), s(L(u'^)), and AR{L{u'^)) needs only 0{n) operations 
[51 [7] ■ At each PCG/ PBiCGstab iteration, we need to calculate the matrix- 
vector product A^j^kW and A'^^w and solve the system Xuy = b. The vector 
multiplication D~^/'^w can be computed in Oin) operations since is a 

diagonal matrix. L{u^)w can be done in 0{n) operations. For H E C oi 
H E ATZ, Hw, H*Hw, and H'Hw can be calculated in O(nlogn) operations 
by few FCTs or FSTs plus lower order of computations. The system Xoy = b 
can also be solved in 0(ri log n) operations. Therefore, the total cost of each 
PCG/ PBiCGstab iteration is bounded by O(nlogn). 



3.2 Two-dimensional Problems 

We can extend the results in Subsection l3.1l to two-dimensional image deblurring 
problems with different BCs. In the two-dimensional case, we assume that the 
PSF is quadrantally symmetric and normalized. When one imposes the reflective 
BCs on H and the zero Neumann BCs on L{u^), the blurring matrix _ff is a 
block Toeplitz-plus-Hankel matrix with Toeplitz-plus-Hankel blocks, which can 
be diagonalized by the two-dimensional FCTs in O(n^logn) operations [20] . 
For an n^-hy-n? matrix A in the form of 



A = 



/ ^1.1 

^2.1 



1,2 
2,2 



^1," \ 
A 

n. 



(22) 
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where Aij are n-by-n matrices, as defined in [5], the Level-1 cosine transform 
preconditioner ci{A) is given by 



ciiA) 



( c(Ai,i) 

C(^2,l) 



C(^l,2) 
C(^2,2) 

C(A„,2) 



C(^l,") \ 
C(A2,„) 

c(^n,n) / 



and then the Level-2 cosine transform preconditioner is C2(^) = Qc\{Ci^ c\(^A)Q)Q^ , 
where Q be the permutation matrix which satisfies [Q"^A<9]ij;fe,; = [^]fe,i;i.j for 
1 ^ j ^ and 1 < fc, Z < n, i.e., the (i, j)th entry of the (fc, Z)th block of A is 
permuted to the (fc, Z)th entry of the (i, j)th block. 

For the two-dimensional linear equation ([5]), using ci(li) = i7, we define the 
optimal reflective BCs preconditioner for A^k in ([5]) by 



R = H*H + ac2{L{u'')). 



(23) 



For eliminating the possibility of large variations in the coefficient of the 
elliptic operator in ([5]) , we employ the same strategy as in Section 13.11 by the 
diagonal scaling in (fT4|) . Therefore, the scaled reflective BCs preconditioner is 
given by 

Dr^D^-RD^. (24) 

where R is defined in (P^ . Similarly, for the scaled system in ([T5|) . the re- 
flective BCs preconditioner is given by Rd = H*H + ac2{L{u^)), where H = 
Hc2{D~^/^). Let Kh, and denote the eigenvalue matrices oiH, 02(0"^^^) 
and C2{L{u'')), respectively. The preconditioner Rj) in ([Ml) can be written as 

Rd = (C„ «) Cnf{A*HAHA*DJ^D + a A^)(C„ C„), 

and hence it is easily inverted by employing few FCTs in 0{n^ log n) operations. 

Next, we assume the anti-reflective BCs for H and the Neumann BCs or 
anti-reflective BCs for L{u''). Then, we construct the anti-reflective BCs pre- 
conditioners for ([5]). For an n^-by-n^ matrix A in (j22l) . the Level-1 sine-based 
transform preconditioner si{A) is given by 



Sl 



(A) = 



s(^2,i) 



^(^1,2) 
•5(^2,2) 



V s(A„,i) S(A„,2) 



S(A2,„) 
s{An^n) / 



By using the same proof technique of Theorem 3.3 in |20] . we can easily show 
that the Level-2 sine-based transform preconditioner S2{A) is given by S2{A) = 
Qsi{Q^ si{A)Q)Q^ . Notice that the matrix H is the anti-reflective BCs matrix. 
Then, we design the sine-based transform preconditioner for ([5]) by 



M = S2{H)*S2{H) + as2(i(u^')). 
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By employing the diagonal scaling in we define the scaled anti-reflective 

BCs preconditioner Dm = D^^'^MD^/^ for ([5]) and the anti-reflective BCs pre- 
conditioner Md = H*H + ac2{L{u^)), where H = S2(-H")s2(D"^/^), for the 
two-dimensional system (IT5|) . Let Kh, and denote the eigenvalue matri- 
ces of S2{D~^/^) and S2{L{u'')), respectively. Then, the preconditioner Mo 
takes the form 

Md = {Sn (E) Sn){A*HAHA*D^D +aAi){Sn (E) Sn), 

which is computationally attractive via FSTs since any matrix operation can be 
done within 0(rt^ log n) operations. 

Finally, we assume the anti-reflective BCs for H and the Neumann BCs 
for L{u''). For an n^-by-n^ matrix A defined in ([2^ . the Level- 1 reblurring 
preconditioner ARi (A) is given by 



ARi{A) = 



/ AR{Ai^i) Ai?(Ai,2) ••• Ai?(Ai,„) \ 
AR{A2,i) ARiA2,2) ■■■ AR{A2,n) 



V Ai?(A„,i) AR{An,2) ■■■ AR{An,n) ) 



Using the same proof as in [20l Theorem 3.3] , we can easily show that the Level-2 
reblurring preconditioner yli?2(^) is given by Ai?2 ( A) = Q^i?i(Q'^Ai?i(A)Q)Q^. 
Now, we design the reblurring preconditioner AR2{A' ^) for the linear equation 
(|20|). Since H is the anti-reflective BCs matrix, we define a reblurring precon- 
ditioner for the linear equation (I^Ul) as P = H'H + a AR2{L{u^)). Also, the 
reblurring preconditioner with diagonal scaling is given by Dp = D^/^PD^^^ 
and the reblurring preconditioner for the two-dimensional scaled linear system 
(HI]) is Pr, = H'H + aAR2{L{u'')), where H = H ■ AR2{D-^/^). Let Ah, Ad, 
and A^ denote the eigenvalue matrices of H, AR2{D^^/^), and AR2{L{u'')), 
respectively. Then, the two-dimensional preconditioner Pd can be written as 

Pd = {Tn®Tn){A*HAHA*DAD+aAi){Tn®Tny\ 

Again, these two-dimensional preconditioner shows interesting computational 
features since the associated linear systems can be solved within O(n^logn) 
operations. 



4 Asymptotic spectral analysis of the precondi- 
tioned sequences 

In order to study the effectiveness of the proposed preconditioners, we need 
the clustering analysis of the spectrum. Also, localization of eigenvalues is of 
interest when solving (|15p via PCG or ((2T|) by PBiCGstab '1 . Here is a useful 
definition |25| for sequences of matrices {A„} where An has size d„, n positive 
integer, and dk > dq li k > q. 
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Definition 3 A matrix sequence {An} is said to be distributed (in the sense of 
the eigenvalues) as the pair (9,G), or have the distribution function 9, if, for 
any F G Co{C), the following limit relation holds 

hm iVF(A,(A„))-^ / t = {ti,...,td), (25) 

where {Aj(A„)}"^]^ denote the eigenvalues of An and fi{-) is the standard Lebesgue 
measure. In that case we write {A„} {9,G). 

An interesting consequence of the equation (|25|) is that {An} ^\ {d,G) 
imphes that most of the eigenvalues are contained within any e-neighborhood 
of the essential range of 6. That is, the range of 6* is a cluster for the spectrum 
of the sequence {A„}. 

The main observation is that all the matrices considered so far are low- 
rank perturbations of Toeplitz matrices or can be viewed as extracted from 
Generalized Locally Toeplitz (GLT) sequences (see [25] and references therein 
and the seminal work [37]). This observation is very important since every GLT 
sequence has a symbol and this symbol is the spectral distribution function 
in the sense of the latter definition. Furthermore, the class of GLT sequences 
is an algebra of matrix-sequences. Hence, when making linear combinations, 
products, or inverses (the latter operation only when the symbol does not vanish 
on sets of positive measure), the result is a new GLT sequence whose symbol can 
be obtained via the same operations on the original symbols, as those performed 
on the matrices. Therefore, a particular case is that of the preconditioned 
matrices can be seen again as extracted from a GLT sequence whose symbol is 
the ratio of the symbols: here the numerator is the symbol of the original matrix 
sequence and the denominator is the symbol of the preconditioning sequence. 

In this section, according to Definition [3] and since we are interested in 
asymptotic estimates, we are forced to indicate explicitly the parameter n which 
uniquely defines the size of the associated matrix. First, we discuss in detail the 
case of reflective BCs. When considering Bn = L{u''), it is well-known [351 that 

{Bn} {a{x)w{t),G) , G = 1} X [0,27r]^ 

d 

a(x) ^ — L w(t) ^ y2 (2-2 cos(t^)). 

On the other hand, c(i3„) ~a (aw(t),G), where a is a constant and in fact 
it is the mean of the function a{x): a — ^(q) /n cl{x) dx . The sequence 
{c{Bn)^^Bn} is clustered at one only if the sequence {Bn — c(i?„)} is clus- 
tered at zero. Since {i?„ — c(i3„)} ^\ {{a{x) — a)w{t),G), the optimal co- 
sine preconditioner is effective only if the function a{x) has no large varia- 
tion. To obtain a clustering preconditioner, a diagonal scaling has to be in- 
troduced. Indeed, the preconditioner diag(i?„)^/^c(i3„)diag(B„)"'^/^ is such that 
{diag(i3„)^/^c(i?„)diag(i?„)^/^} {a{x)w{t),G) due to the algebra stucture 
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of GLT sequences, and hence the preconditioned sequence is clustered at one. 
In our case, the coefficient matrix 

A,, = H*H + aL{u'') 

is the sum of an integral approximate operator and an approximate elhptic 
differential operator. We note that {An} {\h{t)\'^ + aa{x)w{t),G), where 
h is the symbol of the PSF defined in © for the ID case and similarly can 
be defined for d > 1 (the entries of the PSF are the Fourier coefficients of 
h). An effective preconditioner has to consider both terms which consitute the 
matrix An- This is the aim of the preconditioner i?„ defined in (|13l) and (|23p 
for the ID and 2D case, respectively. We have {i?„} ^\ {\h{t)\'^ + aaw{t),G) 
and so {An — Rn} ~a ((a(a;) — a)aw{t),G). In this case, we can not apply 
a diagonal scaling to c{L{u^)) because otherwise we loose the computational 
efficiency, the matrix H*H + adiag(L(M''))^/^c(L(w''))diag(L(w''))^/^ can not 
be diagonalized by discrete cosine transforms. Therefore, we have to apply to 
Rn a diagonal scaling which should be act like diag(L(u'^))^/^ on c{L{u^)), while 
it should be no affect the term H*H . Unfortunately, since we have a diagonal 
scaling we can not apply the diagonal scaling only to a term of the sum. To 
balance the contribution of the two terms, the diagonal scaling is defined by the 
matrix D in ((Ti)) which leads to the preconditioner Da- We have {(£)/{)„} ~;)^ 
((1 + aa{x)){\h{t)\^ + aw{t)),G) and hence the preconditioned sequence is not 
clustered at one, even if for values of a used in the considered applications it 
shows an optimal behaviour (see Figure [5]) . We recall that the clustering is a 
useful property but it is not strictly necessary for the optimality of the related 
preconditioned Krylov method: for instance in the Hermitian positive definite 
case and when dealing with the PCG iterations, the spectral equivalence is 
sufficient. Since An is similar to R~^An, with An = D~^/^AnD~^/'^, the 
use of the preconditioner to the linear system ([5]) is equivalent to apply the 
preconditioner R to the scaled linear system (fTSi) . However, the scaling of the 
linear suggest to use a cosine preconditioner different from R that is Ru which 
is more effective for large values of a (see numerical results in Section [5]). 

Remark 4 For small values of a, i.e., when few regularization is required, the 
three preconditioners R, Dj^ and Rd have a .similar behaviour. Moreover, when 
a goes to zero the effectiveness of the proposed preconditioners increases because 
the preconditioners and the original coefficient matrix An all tend to H* H . 

For concluding this section, we note that in the case of anti-reflective BCs 
similar considerations can be done. The main difference is when we consider the 
reblurring strategy. However, using the results in |14] . the nonsymmetric case 
can be considered as well since the antisymmetric part has trace norm (sum of 
all singular values) bounded by a pure constant independent of n. Therefore 
the spectral distribution is governed by the symmetric part which is dominant 
as discussed in Section 3.3 of [3]. 
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(a) True signal 



(b) Observed signal 



Figure 1: True and observed signals 



5 Numerical Tests 

We will solve the problem by the fixed point method ([S]) with the operator 
T-L approximated by using different BCs and with the matrix C imposed by zero 
Neumann BCs or anti-reflective BCs. The algorithm was implemented in MATLAB 
7 . 10 and run on a PC Intel Pentium IV of 3.00 GHZ CPU. We shall show the 
effectiveness of the proposed preconditioners for the signal/image deblurring 
and also give a comparison of the quality of the restored signals/images with 
different BCs. 

In our test, for simplicity, we choose initial guess vP — v for the FP algorithm. 
We shall solve ([5]) by the PCG method when the Neumann BCs imposed on 
L{u^) and solve ([5]) by the PBiCGstab method when the anti-reflective BCs 
imposed on L{u^). Also, we solve (|20|) by the PBiCGstab method. The initial 
guess for the PCG and PBiCGstab methods in kth FP iteration is chosen to 
be the [k — l)th FP iterate. The PCG and PBiCGstab iterations are stopped 
when the residual vector of the linear systems ^ and (PH)) at the fcth iteration 
satisfies ||?'fc||2/||?'o||2 < tol^ where tol is set to 10^^ and 10^^ in the ID and 2D 
case, respectively. 

5.1 ID case: Signal Deblurring 

In our experiments, we suppose the true signal u is given as in Figure [IJa). The 
two vertical lines shown in Figure [TJa) denote the field of view (i.e., [0.1,0.9]) 
of our signal and the signal outside the two vertical lines can be approximated 
by different BCs. The true signal is blurred by the symmetric out of focus PSF: 



where c is the normalization constant such that hi — I and m{n) is the 
center of the PSF which depends on n so that the restored signal lies in the 




(26) 
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Table 1: Average number of PCG/PBiCGstab iterations per FP step varying 
a, with n = 203 and /3 = 0.1. 

interval [0.1,0.9]. A Gaussian noise ij with noise-to-signal ratio |l?7||2/||-ffu|!2 is 
added to the blurred signal. We consider the true signal is blurred by the out 
of focus PSF and then added the Gaussian noise with the noise levels 1%, i.e., 
||ry||2/||i?u||2 = 0.01. Figure mb) show the observed signal. 

We now show that the proposed preconditioners are effective for solving ([5]) 
and (|20p with different BCs. In our numerical experiment, the FP iteration 
is stopped when — u'^^^ Ib/lju'^lb < lO-^. We will concentrate on the 
performance of different choices of preconditioners for various of parameters a, 
j3, and n. 

In Tables [T] and [21 we report the average number of iterations per FP it- 
eration, where iV, / and D denote the number of FP steps, no preconditioner 
and the diagonal scaling preconditioner, respectively. According to Remark IH 
the effectiveness of the proposed preconditioners increases when a decreases. 
Moreover, decreasing a all the proposed preconditioners become equivalents, 
explicitly the PCG/PBiCGstat converges in about the same number of itera- 
tions. 

We note that anti-reflective BCs usually require lesser steps and lesser PCG/ 
BiCGstab iterations per FP step when compared with reflective BCs. This 
shows that the improvement in the model also leads to an improvement in the 
global computational complexity of the numerical methods. This is more evi- 
dent for the optimal restoration since antireflective BCs require a regularization 
parameter a smaller than the reflective BCs (Figures HHS]). 

Figure[2]describes the average PCG/PBiCGstab iterations per FP step vary- 
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Table 2: Average number of PCG/PBiCGstab iterations per FP step varying 
l3, with n = 203 and a = 0.001. 

ing n. We note that the preconditioners with a diagonal scaling show an optimal 
behavior. 

In Figure [3l we present the restored signals for varying /?, e.g., by solving 
system ([5]) when anti-reflective BCs have been imposed. As expected, the re- 
covered signals become shaper when the value of (3 is smaller, the value (3 — 0.1 
gives a restoration sufhciently good anyway. 

We can easily observe from Tables [TJ2] and Figure [5] that the proposed pre- 
conditioners with a diagonal scaling are the most effective preconditioner when 
varying parameters a, /3, and n. Finally, we remark that in all our tests, the pro- 
posed algorithm needs the same number of FP steps for the no-preconditioner 
and preconditioned cases and ||(7(m'^)||2 tends to 0(10-^) or 0(10"^) at the final 
FP iterate. 

To check the quality of restored signals by using different BCs, in Figure 
|4] we show the relative restoration error (RRE), \\ua ~ "lU/llulb, where Ua is 
the total variation regularized solution of the true signal u, versus the regular- 
ization parameter a. Figure [5] gives the restored signals with optimal value of 
the parameter a, where aopt, Re., Fp., and It. denote the optimal value of 
the parameter a, the minimal RRE, the number of FP steps, and the average 
number of PCG/BiCGstab iterations per FP step, respectively. 

From Figure [5] we argue that anti-refiective BCs lead to the most accurate 
restored signals with less significant ringing effects at the edges and less PCG 
iterations per FP step, when compared with reflective BCs. Moreover, thanks 
to the improvement in the model of the problem, antireflective BCs require 
lesser regularization than reflective BCs. This implies a smaller aopt and hence 
a small number of PCG/BiCGstab iterations per FP step, while the number of 
FP iterations remains about the same. 
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o = I, X = D, + = R, . = D * = R 





(a) 



(b) 



o = I, X = D, + = P. 





(c) 



Figure 2: Average number of PCG/PBiCGstab iterations per FP step for various 
n, with a = 10~^ and /3 = 0.1. (a) Reflective BCs. (b) Anti- Reflective BCs 
with sine preconditioner. (c) Anti-Reflective BCs with reblurring by imposing 
zero Neumann BCs on C (d) Anti-Reflective BCs with reblurring by imposing 
Anti-Reflective BCs on C. 
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Figure 3: Restorations for Anti-Reflective BCs based ([5]) with n = 203, a — 10' 
varying /3. 
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Figure 4: The RRE versus the regularization parameter a for different BCs. 
5.2 2D case: Image Deblurring 

In this section, we apply the proposed preconditioners to image restoration "with 
different BCs. Suppose the true images are bhirred by the Gaussian blur and 
then suppose that a "white Gaussian noise rj "with the noise levels 0.1% is added. 
Figure [5] sho'ws the true and observed images. 

In our numerical tests, the FP iteration is stopped "when — m'''^"'^||2/||u'''||2 
< 10~^ and the maximal number of FP steps is set to be 100. We fix /? = 0.01 
and only focus on the performance of different choices of preconditioners for 
varying a. 

In Table [3] the number of iterations is displayed for solving ([5|) and ([20|) "with 
different BCs and various values of a, "where N and / mean the number of FP 
iterations and no preconditioner, respectively. Here, we only give the average 
number of CG/BiCGstab iterations per FP step. Table [3] suggests that the 
preconditioners Xjj, "with X £ {i?, M, P}, are very effective matrix approxi- 
mations for all values of a, -while the preconditioners Dx do not "work -well for 
large values of a especially for antireflective BCs. Ho^wever, for antirefiective 
BCs a good choice for a is in the interval [10~'^,10~^] and in this case both 
choices Xd and Dx have a similar behaviour. We note that the number of FP 
iterations decrease "with a, so if "we have a good model that requires a lo"wer 
regularization we obtain a gain also in terms of the computational cost of the 
■whole restoration procedure. In all our tests, it is sho"wn that ||5(u'°)||2 tends to 
0(10-3) or 0(10"'') at the final FP iterate. 

Next, we shall check the quality of restored images by using different BCs. 
Figure [7| describes the relative restoration error (RRE) \\ua — u||2/||ti||2 versus 
the regularization parameter a. Figure [8] presents the restored images "with 
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I = 0.047149, Re. = 0.2112. Fp. = 28, It. = 41 0^^^ = 0.00233, Re. = 0.070566, Fp. = 28, It. = 25 
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gure 5: Restored signals with different BCs. Here n = 203 and /3 = 0.1 



True image Observed image 




Figure 6: True and observed images 
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Table 3: Average number of CG/BiCGstab iterations per FP step for varying 
a. Here, (3 = 0.01 (* means that the method does not converge). 
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Figure 7: The RRE versus the regularization parameter a for different BCs. 
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optimal value of the parameter a. Like in the ID case, Figure |5] shows that the 
anti-reflective BCs lead to better restored images and shaper edges with a lower 
computational cost than the reflective BCs (note the high reduction in the FP 
iterations being aopt smaller). 

6 Conclusions 

In this paper, we have considered the effect of reflective and anti-reflective BCs, 
when regularizing blurred and noisy images via the total variation approach. 
In particular, we have studied some preconditioning strategies for the linear 
systems arising from the FP iteration given in . In the case of anti-reflective 
BCs we have also considered a comparison with the reblurring idea proposed in 
[12] • We recall that the latter has been shown effective when combined with the 
Tikhonov regularization and here one of the issues was to verify that reblurring 
and total variation can be combined satisfactorily. Furthermore, the optimal 
behavior of our preconditioners has been validated numerically. 

Beside the computational features of the preconditioning techniques, we 
stress the improvement obtained via anti-reflective BCs both in terms of re- 
construction quality and reduction of the computational cost. In fact, the pre- 
cision of such BCs was already known in the relevant literature (see PH,[TU1[^ 
and references there reported). However, this is the first time that the anti- 
reflective BCs have been combined with a sophisticated regularization method, 
where the use of fast transforms is very welcome for saving computational cost. 
The precision of the antireflective BCs implies also a further reduction of the 
computational cost over the other BCs (see numerical results in Section [5]) be- 
cause they require a smaller regularization parameter a and the effectiveness of 
the proposed preconditioners increases reducing a. 

Potential lines of interest for future investigations could include the use of 
anti-reflective BCs in the promising split Bregman method recently proposed in 
[13] and a more precise clustering analysis of the preconditioning sequences, in 
the spirit of Section S) 
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